Initiating

library(oro.dicom)
library(tidyverse)
library(plotly)
library(data.table)
library(htmlwidgets)
library(coreCT)
library(magick)
library(raster)
library(doParallel)
library(oro.dicom)
library(magick)
library(imager)
ref<-1000*data.frame(hmax=(4.64 +0.77),hmin=(4.64 -0.77), mmax=(3.14 + .65), mmin =(3.14- .65 ))

osis<-path.expand("E:/projects/kaggle/osic-pulmonary-fib/osic-pulmonary-fibrosis-progression")
osistrain<-path.expand("E:/projects/kaggle/osic-pulmonary-fib/osic-pulmonary-fibrosis-progression/train")
osistest<-path.expand("E:/projects/kaggle/osic-pulmonary-fib/osic-pulmonary-fibrosis-progression/test")
osisimg<-path.expand("E:/projects/kaggle/osic-pulmonary-fib/img")

setwd(osis)

#invalids<-"ID00011637202177653955184"   #https://www.kaggle.com/rashmibanthia/corrupt-files

trd<-read.csv("train.csv")  %>% dplyr::rename(Week=Weeks)
#subset(Patient!= invalids)
ted<-read.csv("test.csv") %>% dplyr::rename(Week=Weeks)
invalids<-c("ID00011637202177653955184","ID00052637202186188008618", "ID00128637202219474716089","ID00132637202222178761324")

dd<-list.files(osistrain)
dd <- setdiff(dd,invalids)

*** load dicom img and convert to HU units

x=dd[1]
There were 32 warnings (use warnings() to see them)
read_x<-function(x,thr){
  path<-file.path(osistrain,x)
    s<-readDICOM(path)
    #img<-s$img[[index]] 
    
    ct.slope <- unique(extractHeader(s$hdr, "RescaleSlope"))
    ct.int   <- unique(extractHeader(s$hdr, "RescaleIntercept")) 
    Ct.n<- extractHeader(s$hdr, "InstanceNumber")
    ss <- lapply(s$img, function(x) x*ct.slope + ct.int) 
   # ff<-function(x){x [ x < -2000 ]<- NA;return(x)}
    # trim background
    ss<-  lapply(ss,function(x){x[x< thr]<-F;return(x)})
    #reorder:
    fo<-function(i){
      index<-which(Ct.n ==i)
      return(ss[[index]])
    }
    so<-lapply(1:length(ss),fo)
    return(so)
}
so<-read_x(x,-2000)
    

**imager

l.im<-lapply(t(so[1:3]),as.cimg)
img<-as.imlist(l.im) ; names(img)<-NULL
plot(img,axes = F)

img<-t(so[[7]])



im<-as.cimg( img )
plot(im)

*** more imager processing

px<-(im < - 600 & im > -860 )

plot(px)

plot(im)
highlight(px)

NA
NA
plot(im)
px.flood(im,78,296,sigma=.21) %>% highlight



pal <- colorRamp(c("black","white"))

fig <- plot_ly(z = ~ so[[7]], colors = pal) %>%
add_surface(
  contours = list(
    z = list(
      show=TRUE,
      usecolormap=TRUE,
      highlightcolor="#ff0000",
      project=list(z=TRUE)
    )
  )
)
fig <- fig %>% layout(
  scene = list(
    camera=list(
      eye = list(x=1.87, y=0.88, z=-0.64)
    )
  )
)

fig

NA

** new threshould to focus on lungs

sfail<-read_x(x,-600)
'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2
pal <- colorRamp(c("black","white"))

fig <- plot_ly(z = ~ sfail[[7]], colors = pal) %>%
add_surface(
  contours = list(
    z = list(
      show=TRUE,
      usecolormap=TRUE,
      highlightcolor="#ff0000",
      project=list(z=TRUE)
    )
  )
)
fig <- fig %>% layout(
  scene = list(
    camera=list(
      eye = list(x=1.87, y=0.88, z=-0.64)
    )
  )
)

fig

It has cut most of the lungs…that is below -800

**** lets start lung segmentation linear model:

d <- as.data.frame(im)
##Subsamble, fit a linear model
m <- sample_n(d,1e4) %>% lm(value ~ x*y,data=.) 
##Correct by removing the trend
im.c <- im-predict(m,d)
out <- threshold(im.c)
plot(out)

 out<-threshold(im,"auto")
out <- clean(out,3) %>% imager::fill(7)
plot(im)
highlight(out)

NA
NA
im[ out]<- NA
plot(im)

 plot_ly(
            z = so[[7]],
            #colorscale = list(c(0,0.5,1),c("blue", "white", "red")),
            colors=pal,
            type = "heatmapgl"
            )

NA
NA
a<-kmeans(out,2)
colors<- a$centers[a$cluster[1,],]
im<-as.cimg(so[[7]]) %>% imrotate(-90)


detect.edges <- function(im,sigma=1)
    {
        isoblur(im,sigma) %>% imgradient("xy") %>% enorm %>% imsplit("c") %>% add
    }

edges <- detect.edges(im,2) %>% sqrt 
plot(edges)

pmap <- 1/(1+edges) #Priority inv. proportional to gradient magnitude
plot(pmap,main="Priority map") #Nice metal plate effect! 

NA
NA
seeds <- imfill(dim=dim(pmap)) #Empty image
seeds[170,233,1,1] <- 1 #Background pixel 
seeds[36,255,1,1] <- 2 #Foreground pixel

wt1 <- watershed(seeds,pmap)
plot(wt1,main="Watershed segmentation")

the other lung

seeds <- imfill(dim=dim(pmap)) #Empty image
seeds[355,247,1,1] <- 1 #Background pixel 
seeds[56,252,1,1] <- 2 #Foreground pixel

wt2 <- watershed(seeds,pmap)
plot(wt2,main="Watershed segmentation")

NA
NA

now, both


wt <- wt1 + wt2

plot(wt)

NA
NA
wt[wt==3]<-T
wt[wt==4]<-F
lungs<-wt * im
plot(lungs)

NA
NA

dx <- imgradient(im,"x")
dy <- imgradient(im,"y")
grad.mag <- sqrt(dx^2+dy^2)
plot(grad.mag,main="Gradient magnitude")

** lets try kmeans


a<-kmeans(im)
plot(im) %>% highlight(a$centers)
LS0tDQp0aXRsZTogIlIgTm90ZWJvb2siDQpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sNCi0tLQ0KDQpJbml0aWF0aW5nDQoNCmBgYHtyIG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0V9DQpsaWJyYXJ5KG9yby5kaWNvbSkNCmxpYnJhcnkodGlkeXZlcnNlKQ0KbGlicmFyeShwbG90bHkpDQpsaWJyYXJ5KGRhdGEudGFibGUpDQpsaWJyYXJ5KGh0bWx3aWRnZXRzKQ0KbGlicmFyeShjb3JlQ1QpDQpsaWJyYXJ5KG1hZ2ljaykNCmxpYnJhcnkocmFzdGVyKQ0KbGlicmFyeShkb1BhcmFsbGVsKQ0KbGlicmFyeShvcm8uZGljb20pDQpsaWJyYXJ5KG1hZ2ljaykNCmxpYnJhcnkoaW1hZ2VyKQ0KcmVmPC0xMDAwKmRhdGEuZnJhbWUoaG1heD0oNC42NCArMC43NyksaG1pbj0oNC42NCAtMC43NyksIG1tYXg9KDMuMTQgKyAuNjUpLCBtbWluID0oMy4xNC0gLjY1ICkpDQoNCm9zaXM8LXBhdGguZXhwYW5kKCJFOi9wcm9qZWN0cy9rYWdnbGUvb3NpYy1wdWxtb25hcnktZmliL29zaWMtcHVsbW9uYXJ5LWZpYnJvc2lzLXByb2dyZXNzaW9uIikNCm9zaXN0cmFpbjwtcGF0aC5leHBhbmQoIkU6L3Byb2plY3RzL2thZ2dsZS9vc2ljLXB1bG1vbmFyeS1maWIvb3NpYy1wdWxtb25hcnktZmlicm9zaXMtcHJvZ3Jlc3Npb24vdHJhaW4iKQ0Kb3Npc3Rlc3Q8LXBhdGguZXhwYW5kKCJFOi9wcm9qZWN0cy9rYWdnbGUvb3NpYy1wdWxtb25hcnktZmliL29zaWMtcHVsbW9uYXJ5LWZpYnJvc2lzLXByb2dyZXNzaW9uL3Rlc3QiKQ0Kb3Npc2ltZzwtcGF0aC5leHBhbmQoIkU6L3Byb2plY3RzL2thZ2dsZS9vc2ljLXB1bG1vbmFyeS1maWIvaW1nIikNCg0Kc2V0d2Qob3NpcykNCg0KI2ludmFsaWRzPC0iSUQwMDAxMTYzNzIwMjE3NzY1Mzk1NTE4NCIgICAjaHR0cHM6Ly93d3cua2FnZ2xlLmNvbS9yYXNobWliYW50aGlhL2NvcnJ1cHQtZmlsZXMNCg0KdHJkPC1yZWFkLmNzdigidHJhaW4uY3N2IikgICU+JSBkcGx5cjo6cmVuYW1lKFdlZWs9V2Vla3MpDQojc3Vic2V0KFBhdGllbnQhPSBpbnZhbGlkcykNCnRlZDwtcmVhZC5jc3YoInRlc3QuY3N2IikgJT4lIGRwbHlyOjpyZW5hbWUoV2Vlaz1XZWVrcykNCmludmFsaWRzPC1jKCJJRDAwMDExNjM3MjAyMTc3NjUzOTU1MTg0IiwiSUQwMDA1MjYzNzIwMjE4NjE4ODAwODYxOCIsICJJRDAwMTI4NjM3MjAyMjE5NDc0NzE2MDg5IiwiSUQwMDEzMjYzNzIwMjIyMjE3ODc2MTMyNCIpDQoNCmRkPC1saXN0LmZpbGVzKG9zaXN0cmFpbikNCmRkIDwtIHNldGRpZmYoZGQsaW52YWxpZHMpDQpgYGANCg0KDQoNCioqKiBsb2FkIGRpY29tIGltZyBhbmQgY29udmVydCB0byBIVSB1bml0cw0KDQpgYGB7ciBtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPUZBTFNFfQ0KeD1kZFsxXQ0KDQpyZWFkX3g8LWZ1bmN0aW9uKHgsdGhyKXsNCiAgcGF0aDwtZmlsZS5wYXRoKG9zaXN0cmFpbix4KQ0KICAgIHM8LXJlYWRESUNPTShwYXRoKQ0KICAgICNpbWc8LXMkaW1nW1tpbmRleF1dIA0KICAgIA0KICAgIGN0LnNsb3BlIDwtIHVuaXF1ZShleHRyYWN0SGVhZGVyKHMkaGRyLCAiUmVzY2FsZVNsb3BlIikpDQogICAgY3QuaW50ICAgPC0gdW5pcXVlKGV4dHJhY3RIZWFkZXIocyRoZHIsICJSZXNjYWxlSW50ZXJjZXB0IikpIA0KICAgIEN0Lm48LSBleHRyYWN0SGVhZGVyKHMkaGRyLCAiSW5zdGFuY2VOdW1iZXIiKQ0KICAgIHNzIDwtIGxhcHBseShzJGltZywgZnVuY3Rpb24oeCkgeCpjdC5zbG9wZSArIGN0LmludCkgDQogICAjIGZmPC1mdW5jdGlvbih4KXt4IFsgeCA8IC0yMDAwIF08LSBOQTtyZXR1cm4oeCl9DQogICAgIyB0cmltIGJhY2tncm91bmQNCiAgICBzczwtICBsYXBwbHkoc3MsZnVuY3Rpb24oeCl7eFt4PCB0aHJdPC1GO3JldHVybih4KX0pDQogICAgI3Jlb3JkZXI6DQogICAgZm88LWZ1bmN0aW9uKGkpew0KICAgICAgaW5kZXg8LXdoaWNoKEN0Lm4gPT1pKQ0KICAgICAgcmV0dXJuKHNzW1tpbmRleF1dKQ0KICAgIH0NCiAgICBzbzwtbGFwcGx5KDE6bGVuZ3RoKHNzKSxmbykNCiAgICByZXR1cm4oc28pDQp9DQpzbzwtcmVhZF94KHgsLTIwMDApDQogICAgDQpgYGANCioqaW1hZ2VyDQoNCmBgYHtyfQ0KbC5pbTwtbGFwcGx5KHQoc29bMTozXSksYXMuY2ltZykNCmltZzwtYXMuaW1saXN0KGwuaW0pIDsgbmFtZXMoaW1nKTwtTlVMTA0KcGxvdChpbWcsYXhlcyA9IEYpDQoNCmBgYA0KDQoNCmBgYHtyfQ0KaW1nPC10KHNvW1s3XV0pDQoNCg0KDQppbTwtYXMuY2ltZyggaW1nICkNCnBsb3QoaW0pDQpgYGANCg0KKioqIG1vcmUgaW1hZ2VyIHByb2Nlc3NpbmcNCg0KDQpgYGB7cn0NCnB4PC0oaW0gPCAtIDYwMCAmIGltID4gLTg2MCApDQoNCnBsb3QocHgpDQpgYGANCg0KDQoNCg0KYGBge3J9DQpwbG90KGltKQ0KaGlnaGxpZ2h0KHB4KQ0KDQoNCmBgYA0KYGBge3J9DQpwbG90KGltKQ0KcHguZmxvb2QoaW0sNzgsMjk2LHNpZ21hPS4yMSkgJT4lIGhpZ2hsaWdodA0KDQpgYGANCg0KDQoNCg0KDQoNCg0KYGBge3J9DQoNCg0KcGFsIDwtIGNvbG9yUmFtcChjKCJibGFjayIsIndoaXRlIikpDQoNCmZpZyA8LSBwbG90X2x5KHogPSB+IHNvW1s3XV0sIGNvbG9ycyA9IHBhbCkgJT4lDQphZGRfc3VyZmFjZSgNCiAgY29udG91cnMgPSBsaXN0KA0KICAgIHogPSBsaXN0KA0KICAgICAgc2hvdz1UUlVFLA0KICAgICAgdXNlY29sb3JtYXA9VFJVRSwNCiAgICAgIGhpZ2hsaWdodGNvbG9yPSIjZmYwMDAwIiwNCiAgICAgIHByb2plY3Q9bGlzdCh6PVRSVUUpDQogICAgKQ0KICApDQopDQpmaWcgPC0gZmlnICU+JSBsYXlvdXQoDQogIHNjZW5lID0gbGlzdCgNCiAgICBjYW1lcmE9bGlzdCgNCiAgICAgIGV5ZSA9IGxpc3QoeD0xLjg3LCB5PTAuODgsIHo9LTAuNjQpDQogICAgKQ0KICApDQopDQoNCmZpZw0KDQpgYGANCg0KKiogbmV3IHRocmVzaG91bGQgdG8gZm9jdXMgb24gbHVuZ3MNCg0KYGBge3J9DQpzZmFpbDwtcmVhZF94KHgsLTYwMCkNCnBhbCA8LSBjb2xvclJhbXAoYygiYmxhY2siLCJ3aGl0ZSIpKQ0KDQpmaWcgPC0gcGxvdF9seSh6ID0gfiBzZmFpbFtbN11dLCBjb2xvcnMgPSBwYWwpICU+JQ0KYWRkX3N1cmZhY2UoDQogIGNvbnRvdXJzID0gbGlzdCgNCiAgICB6ID0gbGlzdCgNCiAgICAgIHNob3c9VFJVRSwNCiAgICAgIHVzZWNvbG9ybWFwPVRSVUUsDQogICAgICBoaWdobGlnaHRjb2xvcj0iI2ZmMDAwMCIsDQogICAgICBwcm9qZWN0PWxpc3Qoej1UUlVFKQ0KICAgICkNCiAgKQ0KKQ0KZmlnIDwtIGZpZyAlPiUgbGF5b3V0KA0KICBzY2VuZSA9IGxpc3QoDQogICAgY2FtZXJhPWxpc3QoDQogICAgICBleWUgPSBsaXN0KHg9MS44NywgeT0wLjg4LCB6PS0wLjY0KQ0KICAgICkNCiAgKQ0KKQ0KDQpmaWcNCmBgYA0KDQpJdCBoYXMgY3V0IG1vc3Qgb2YgdGhlIGx1bmdzLi4udGhhdCBpcyBiZWxvdyAtODAwDQoNCg0KDQogKioqKiBsZXRzIHN0YXJ0IGx1bmcgc2VnbWVudGF0aW9uDQpsaW5lYXIgbW9kZWw6DQoNCmBgYHtyfQ0KZCA8LSBhcy5kYXRhLmZyYW1lKGltKQ0KIyNTdWJzYW1ibGUsIGZpdCBhIGxpbmVhciBtb2RlbA0KbSA8LSBzYW1wbGVfbihkLDFlNCkgJT4lIGxtKHZhbHVlIH4geCp5LGRhdGE9LikgDQojI0NvcnJlY3QgYnkgcmVtb3ZpbmcgdGhlIHRyZW5kDQppbS5jIDwtIGltLXByZWRpY3QobSxkKQ0Kb3V0IDwtIHRocmVzaG9sZChpbS5jKQ0KcGxvdChvdXQpDQoNCmBgYA0KDQoNCg0KYGBge3J9DQogb3V0PC10aHJlc2hvbGQoaW0sImF1dG8iKQ0Kb3V0IDwtIGNsZWFuKG91dCw1KSAlPiUgaW1hZ2VyOjpmaWxsKDcpDQpwbG90KGltKQ0KaGlnaGxpZ2h0KG91dCkNCg0KDQpgYGANCg0KYGBge3J9DQppbW88LWltDQppbW9bIG91dF08LSBOQQ0KcGxvdChpbW8pDQpgYGANCg0KYGBge3J9DQogcGxvdF9seSgNCiAgICAgICAgICAgIHogPSBzb1tbN11dLA0KICAgICAgICAgICAgI2NvbG9yc2NhbGUgPSBsaXN0KGMoMCwwLjUsMSksYygiYmx1ZSIsICJ3aGl0ZSIsICJyZWQiKSksDQogICAgICAgICAgICBjb2xvcnM9cGFsLA0KICAgICAgICAgICAgdHlwZSA9ICJoZWF0bWFwZ2wiDQogICAgICAgICAgICApDQoNCg0KYGBgDQoNCg0KDQpgYGB7cn0NCmE8LWttZWFucyhvdXQsMikNCmNvbG9yczwtIGEkY2VudGVyc1thJGNsdXN0ZXJbMSxdLF0NCg0KDQpgYGANCg0KDQpgYGB7cn0NCmltPC1hcy5jaW1nKHNvW1s3XV0pICU+JSBpbXJvdGF0ZSgtOTApDQoNCg0KZGV0ZWN0LmVkZ2VzIDwtIGZ1bmN0aW9uKGltLHNpZ21hPTEpDQogICAgew0KICAgICAgICBpc29ibHVyKGltLHNpZ21hKSAlPiUgaW1ncmFkaWVudCgieHkiKSAlPiUgZW5vcm0gJT4lIGltc3BsaXQoImMiKSAlPiUgYWRkDQogICAgfQ0KDQplZGdlcyA8LSBkZXRlY3QuZWRnZXMoaW0sMikgJT4lIHNxcnQgDQpwbG90KGVkZ2VzKQ0KcG1hcCA8LSAxLygxK2VkZ2VzKSAjUHJpb3JpdHkgaW52LiBwcm9wb3J0aW9uYWwgdG8gZ3JhZGllbnQgbWFnbml0dWRlDQpwbG90KHBtYXAsbWFpbj0iUHJpb3JpdHkgbWFwIikgI05pY2UgbWV0YWwgcGxhdGUgZWZmZWN0ISANCg0KDQpgYGANCg0KYGBge3J9DQpzZWVkcyA8LSBpbWZpbGwoZGltPWRpbShwbWFwKSkgI0VtcHR5IGltYWdlDQpzZWVkc1sxNzAsMjMzLDEsMV0gPC0gMSAjQmFja2dyb3VuZCBwaXhlbCANCnNlZWRzWzM2LDI1NSwxLDFdIDwtIDIgI0ZvcmVncm91bmQgcGl4ZWwNCg0Kd3QxIDwtIHdhdGVyc2hlZChzZWVkcyxwbWFwKQ0KcGxvdCh3dDEsbWFpbj0iV2F0ZXJzaGVkIHNlZ21lbnRhdGlvbiIpDQoNCg0KDQoNCg0KYGBgDQp0aGUgb3RoZXIgbHVuZw0KDQpgYGB7cn0NCnNlZWRzIDwtIGltZmlsbChkaW09ZGltKHBtYXApKSAjRW1wdHkgaW1hZ2UNCnNlZWRzWzM1NSwyNDcsMSwxXSA8LSAxICNCYWNrZ3JvdW5kIHBpeGVsIA0Kc2VlZHNbNTYsMjUyLDEsMV0gPC0gMiAjRm9yZWdyb3VuZCBwaXhlbA0KDQp3dDIgPC0gd2F0ZXJzaGVkKHNlZWRzLHBtYXApDQpwbG90KHd0MixtYWluPSJXYXRlcnNoZWQgc2VnbWVudGF0aW9uIikNCg0KDQpgYGANCm5vdywgYm90aA0KDQpgYGB7cn0NCg0Kd3QgPC0gd3QxICsgd3QyDQoNCnBsb3Qod3QpDQoNCg0KYGBgDQoNCmBgYHtyfQ0Kd3Rbd3Q9PTNdPC1UDQp3dFt3dD09NF08LUYNCmx1bmdzPC13dCAqIGltDQpwbG90KGx1bmdzKQ0KDQoNCmBgYA0KDQoNCg0KDQoNCg0KDQoNCg0KDQpgYGB7cn0NCg0KZHggPC0gaW1ncmFkaWVudChpbSwieCIpDQpkeSA8LSBpbWdyYWRpZW50KGltLCJ5IikNCmdyYWQubWFnIDwtIHNxcnQoZHheMitkeV4yKQ0KcGxvdChncmFkLm1hZyxtYWluPSJHcmFkaWVudCBtYWduaXR1ZGUiKQ0KDQpgYGANCg0KDQoNCg0KDQoqKiBsZXRzIHRyeSBrbWVhbnMNCg0KDQpgYGB7cn0NCg0KYTwta21lYW5zKGltKQ0KcGxvdChpbSkgJT4lIGhpZ2hsaWdodChhJGNlbnRlcnMpDQoNCg0KYGBgDQpgYGB7cn0NCg0KYGBgDQoNCg==